A/ASA -TW- 2U702 


NASA Technical Memorandum 86708 


NASA-TM-86708 19860020159 


An Outflow Acoustic Boundary 
Condition for Internal Duct Flows 

Marianne Mosher 


February 1986 



I'viAR 1 7 198b 


l-Ai VO - w v . „ : (' • ; G £N T ER 
L! Li IaAI7Y, NASA 

HAMPTON, VIRGINIA 


MSA 

National Aeronautics and 
Space Administration 






NASA Technical Memorandum 86708 


An Outflow Acoustic Boundary 
Condition for Internal Duct Flows 

Marianne Mosher, Ames Research Center, Moffett Field, California 


February 1986 


NASA 

National Aeronautics and 
Space Administration 

Ames Research Center 

Moffett Field, California 94035 


A/<&>-296<3/^ 




SYMBOLS 


c — speed of sound 

Cj = complex harmonic amplitude of the jth duct mode 

dpd = vector of normal derivative of harmonic pressure at 

downstream boundary 

dp u = vector of normal derivative of harmonic pressure at up- 
stream boundary 

e = 2.718282 

f = source term in Helmholtz equation 

i = y/^\ 

j — index for eigenvalues in x direction 

k = wave number, - 

k t = downstream x wave number (eigenvalue) for the jth duct 
mode 

k~ — upstream x wave number (eigenvalue) for the jth duct 
mode 

kd = diagonal matrix of duct downstream eigenvalues 

k u = diagonal matrix of duct upstream eigenvalues 

/ = index for control points 

M = freestream wind tunnel Mach number, Uoo/c 
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normal to the duct surface, positive toward duct center 
indicies for eigenvalues in y and :z direction 
normalizing constants for eigenmodes 
acoustic pressure 

complex pressure coefficient for single frequency, P = '3?(e -,u, *p) 

vector of harmonic pressures on downstream outflow 
boundary 

harmonic pressure amplitude of jth duct eigenmode for 
downstream propagation 

vector of harmonic pressures on upstream outflow 
boundary 

harmonic pressure amplitude of jth eigenmode for up- 
stream propagating wave 

surface of control volume 

downstream outflow control surface 

upstream outflow control surface 

duct wall control surface 

componant of acoustic velocity normal to surface 
axial componant of uniform flow in duct 
control volume in duct 
coordinates 
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acoustic source position 

acoustic impedance of a surface 

boundary condition coefficient 

3.141593 

duct eigenmodes 

matrix of duct eigenvectors 

acoustic frequency 
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AN OUTFLOW ACOUSTIC BOUNDARY CONDITION 


FOR INTERNAL DUCT FLOWS 

Marianne Mosher 
Ames Research Center 


SUMMARY 


A boundary condition for the linear acoustic equation has been developed that 
allows the acoustic pressure waves to propagate out of the computational domain 
boundary, just as they would propagate in an infinitely long duct. The problem 
is divided into two domains: numerical and analytical. The boundary condition 
provides a matching of the two domains. Examples show this method works well 
in an acoustic panel program for a model problem (simple source in a rectangular 
duct with several propagating modes present). This report describes the bound- 
ary condition so that it can be used with various duct geometries and numerical 
methods. 


INTRODUCTION 


When computing acoustic properties in large or infinite regions, it is often de- 
sirable to restrict the computation to a subregion. This requires a condition to 
be imposed at an artificial boundary that will represent the acoustic field in the 
region outside the computational domain. The boundary condition must be very 
accurate for acoustic problems, otherwise the acoustic waves will reflect off the ar- 
tificial boundary back into the computational domain. Acoustic fields in ducts are 
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particularly sensitive; a poor boundary condition can change the result by 100% or 
more. 

Unlike a wall boundary, no simple relation exists between the acoustic pressure 
and its derivative at the artificial outflow boundary. The relationship between the 
two depends on the solution; both on the acoustic source, and how the acoustic 
field propagates in the duct away from the control volume. 

The simplest boundary condition that might be used is a linear relation between p 
and its derivative: dp/dn = ap. This would allow a plane wave with wave number 
in the x-direction of k x = —ia to pass through the boundary without spurious 
reflections. All other wave numbers would be partially reflected. 

In reference 1, Bayliss and Turkel present a series of boundary conditions that 
allow several modes to pass through the boundary. Their boundary condition is 

j=jmax ry 

n ■(^-**i)p = ° (!) 

When all of the infinite terms are included, this exactly models the acoustic waves 
propagating in the duct away from the source. At least one term of the product 
must be included for each propagating wave that exists at the boundary. When 
many terms are used it becomes difficult to apply this boundary condition as it 
involves many high order derivatives. This paper develops a boundary condition 
that involves only one normal derivative. 

The following method may be used for linear acoustic duct problems solved in 
the frequency domain. At the boundary, the acoustic field for one frequency can be 
expressed as a linear combination of eigenmodes, truncated to a finite number of 
propagating eigenmodes and decaying eigenmodes. The eigenvalues and eigenmodes 
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depend on the duct shape, wall boundary, and duct terminations; thus they can be 
determined independently of the source. An equation can be constructed which 
relates p and dp/dn at the boundary. The discrete version of this equation can 
be used as the boundary condition and combined with the discrete equations to be 
solved. Solving the system then produces the amplitudes of the eigenmodes. 

MATHEMATICAL MODEL 


For subsonic flow with Mach number M in the +x direction the equation describ 


ing the acoustic field is the convected Helmholtz equation: 


V 2 p - M 2 ^JJ + 2ikM~ + A; 2 p = /(x, y, z, k) 


■i 


(2) 


where f is the source. The boundary condition on the wall is: 


dp 

ap + -r- = 0, on S 

on 


( 3 ) 
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Figure 1. - Control Volume 
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Figure 1 shows the duct geometry. 


The homogeneous solution consists of acoustic modes propagating downstream 
in the (+x) direction and propagating upstream in the (-x) direction . 


P + = E ct0+( y , *)<>•'* y (■ > x a (4a) 

3=0 

p" = E ( 4fe ) 

3=0 

The functions rpf are the duct modes for the eigenvalues A;* . 


Any linear acoustic field in a constant cross-section duct can be represented by 
a linear combination of duct modes. Far away from all acoustic sources, a finite 
number of modes represents the acoustic field well. The amplitude of most modes 
decays exponentially with distance away from the source. Thus, at the outflow 
boundaries, the sum of a finite number of modes represents the acoustic pressure 
well. For example, at the downstream boundary, where x > x s for all the acoustic 
sources, the harmonic acoustic pressure may be written, 


p(x, !) , 2 )=’f C ,0 y (y, 2 )e'‘. t '-- ) 


Jmax 

E c ;V>;(y, 


*) 


e «*+ (*-*.) 


(5) 


3=1 3=1 

where kj is the jth eigenvalue, ipj is the jth eigenmode, and c, is the amplitude of 
the jth mode. With this restriction to the eigenmodes of the physical duct, both the 
harmonic acoustic pressure and its normal derivative can be expressed in terms of 
the same variables, p dj = e' k > ( z "*“ -I ') c + 5 the complex amplitudes of the j max modes 
at the boundary x = x max . 


Jmax 


P (Xmax, y, Z) » E Vv(y> Z )P dj 

3=1 


(6a) 
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(6 6 ) 


dp { x max>V^ z ) dp(x mal , J/ ? z ) _ 

d7 “ dx * 

Likewise at the upstream boundary with x < x s 


Jmax 

£ 


j=i 


ikftj}j[y,z)i> d . 


3 m ax 

p{x min ,y,z) « 52 
j'=i 

dp (Xmm»y,z) _ dp(j mm? y ? 2) 

dn dx 


0i(y.*)^u y 

3max 

« - 12 

i=i 


(6c) 

(6d) 


Next, consider the discretized variables: evaluate the quantities at discrete (y,z) 
points on the boundary. All the relations in equation (6) can be represented by 
matrices: 


where, 


Pu = ®P U . 

(7a) 

dpu = *®kJi Ui 

(76) 

Pd = 

(7c) 

dp<i = 

(7 d) 

p u = {—p{x m in,yi,Zi)...) T 

(8a) 

^ dp(x m j n , yi, zi) 
dn 

(86) 

Pd = (-P (*ma*.yj,«l)-) T 

(8c) 

dp, = (... 3n ...) 

(id) 


When the number of control points on the outflow plane equals the number of 
modes included in the representation, p can be eliminated by direct inversion to 
obtain a matrix equation relating p and dp. 

i\yk u \P -1 p u — dp = 0 (9a) 
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t^kd® - dp = 0 (96) 

When the number of control points is greater than the number of modes, a least 
squares solution for p gives: 


i¥k a [® T *]” 1 ¥ T p u - dp = 0 

(9c) 

f¥k 11 [<P T ®]" 1 ¥ T p (1 - dp = 0 

(M) 


Outside the control volume a sum of eigenmodes represents the acoustic field. 

2 max 

P = E Vb(y,z)e ,<: ' {z ~ Xmax) p d .,x > Xmax (lOa) 

i = i 

2 max __ 

P = E $i(v> z)e~ ik 7 (x - Xmin) p u .,x < x min (106) 

3=1 

The boundary condition matches the acoustic field in the control volume to an 
analytic solution outside the control volume. 

This form can be used because the acoustic field inside the duct can be expressed 
as a linear sum of modes. Any type of duct can be used as long as the eigenvalues 
and eigenmodes can be determined. Often the duct is simple enough (separable 
coordinates and simple boundary conditions on the wall) that the eigenvalues and 
eigenmodes have analytic expressions. For more complicated ducts, the eigenvalues 
and eigenmodes must be determined numerically. 

All modes with any significant amplitude at the boundary must be included in 
the eigenmode representation. Generally, those will be all propagating modes, and 
possibly a few decaying modes. The boundary should be placed far enough away 
from the sources so that decaying modes will have a chance to decay. 
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EXAMPLES 


Model Problem 


This boundary condition has been incorporated into an acoustic panel method 
and applied successfully to a model problem of a point acoustic source in an infinitely 
long rectangular duct. The duct is bound by y = ±6/2 and z = ±d/2. The 
eigenvalues and eigenmodes may be determined analytically. 

Eigenvalues and Eigenmodes 


By the method of separation of variables as described in reference 2, the solution 
to equation (2) in the rectangular duct with boundary condition in equation (3) 
can be found. Acoustic modes propagating downstream in the +x direction have 
the form: 

n=oo 

m = oo , 

P + = H cs(k lim y)cs(k z „z)e' i > , '"*' ) p+,x > x. 


(11a) 


n — oo 
m=U 


and propagating in the upstream (-x) direction have the form: 

n=oo 

m=oo _ _ 

P _ = J2 cs { k vmy)zs{kz n z)e~' k i ^’’Pj , X < X, 

n=oo 

The function cs(k n ) is defined: 

, . / cos(A; n ), n< 

cs(k n ) = < . ; / 

( sin (k n ),n< 


(116) 


:even 

odd 


( 12 ) 


The boundary conditions give the following trancsendental equations for the eigen- 


values k ym and k zn : 


a = k ym tan(A: t/m -), meven 
a = A: vm tan(^ + k ym ^), modd 


(13a) 

(136) 
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(13c) 


a = k zn ta.n(k zn -), neven 

Jd 

a = k zn tan(^ + nodd 


(13d) 


m7r 

kym — f" £z 


A solution for each eigenvalue will exist for each period of n . The solutions will 
have the form: 

mir 

(14a) 
(146) 

When the wall is close to being a hard wall, a is small, and the solutions give 
approximate values of: 

' m = 0 (15a) 


M7T 

Kzn “j ^yn* 


~ym 


£zm ~ 


b 
8a 
bir 2 m 2 


, m > 0 


2a __ n 

, , n — 0 


(156) 


dit 2 n 2 ’ n > 0. 

The eigenvalues, kj mn and k~ mn for propagation downstream and upstream, are 
found. 

-kM ± \lk* - fl - M 2 U0L.J 2 + (M 2 ) 

— (16a) 


fc + = 


and 


(1 -M 2 ) 


kM ± sfk* - (1 - M 2 )((fc„ m ) 2 + (fc*„) 2 ) 


(T^NPj ' (166) 

The eigenvalues are chosen such that the imaginary part of kj mn > 0 to ensure that 
the decaying waves are propagating away from the source. 


An ordering of the eigenvalues, which includes at least all of the propagating 
modes, maps the m,n indices into the j mode indices. Thus, we have &(l,j) = 
cs(k ym yi)cs(k zn zi) is the mode shape of the yth mode at position (y/, zi), k u (y, j) = kj 
is the jth upstream eigenvalue, and k d{j\j) = kj is the jth downstream eigenvalue. 
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Analytic Solution 


The partial differential equation for a monopole in a duct is: 


v 2 P - m 2 0 + 2ikM % + * 2 p = «(* - *.)% - V.W* ~ *.)• (17) 


The solution will have the following form: 

TO — OO 

n=oo 

P* = £ r - S iKmy) a ( k ™z)fmr,( X ) 


(18) 


m=0 

n=0 


Substituting equation (18) into (17), and solving by a Fourier Transform method 
gives: 

/ » \ / i \ 

. e ±ik£ mn {z-x,)' (jQ) 


r± cs (^ym!/«) cs (^n^ : «) 

J inn 


(l-M‘)N m N n 2(*± m „± &) 

The eigenfunctions are orthogonal with respect to the duct area; however, they are 
not necessarily normal. Let N m and N n be the normalizing constants: 

N„ = /j cs 2 (fc mv y)dy = b + 


2k 


my 


N n = f* d cs 2 {k nz z)dz = d+ (-iy sm ( k nz d ) 


(20a) 

(206) 


Results 


Two cases will be shown for a simple monopole source inside an infinitely long 
rectangular duct without flow. The hard-walled duct has an aspect ratio of 0.7. The 
source is centered longitudinaly in the control volume at x = 0. In the first case, 
the source has a low frequency, producing 1 propagating mode(plane wave) with 
outflow boundaries located at x = ±1.5, and source located off axis at y = 0.09 
and z = 0.08 In the second case, the source has a higher frequency producing 4 
propagating modes with outflow boundaries at x = ±1.2, and source located off 
axis at y = 0.20 and z = 0.15 
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Three solutions are shown for each case: the analytic series solution represented 
by equation (18) (which converges everywhere except in the plane of the source, x — 
0), a computation with the new matched boundary condition, and a computation 
with a simple boundary condition specifying the impedance of the lowest mode at 
the outflow boundary. This simple boundary condition is equivalent to using one 
mode in this matching procedure or one term in the boundary condition of Bayliss 
and Turkel. 

For the the low frequency case, six modes are matched at the outflow boundaries, 
3 eigenvalues in the y direction and 2 eigenvalues in the z direction. Table 1 shows 
the 6 eigenvalues and their coefficients at the outflow boundary. Good agreement 
exists between the analytic coefficients and the coefficients computed with the new 
procedure. The first coefficient computed with the simple boundary condition also 
shows good agreement with the analytic coefficient. Figure 2 shows the real part of 
the pressure solution along 4 longitudinal lines in two horizontal planes. The dotted 
line indicates the location of the outflow boundary. Both computations agree well 
with the analytic solution. Figure 3 shows the amplitude of acoustic pressure. 

For the higher frequency model problem, 20 modes are matched at the outflow 
boundaries, 5 eigenvalues in the y direction and 4 eigenvalues in the z direction. 
Table 2 shows the 20 eigenvalues and their coefficients at the outflow boundary. 
Good agreement exists between the analytic coefficients and those computed with 
the new matched boundary condition procedure, particularly the lower order prop- 
agating modes. The coefficient for the first propagating mode predicted with the 
simple boundary condition fails to match the analytically determined coefficient at 
all. Figures 4 and 5 show the real part and amplitude of the acoustic field respec- 
tively. The computation with the new matched boundary condition agrees well with 
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the analytic solution. The computation with the simple boundary condition differs 
considerably from the analytic solution. 


DISCUSSION 


This method has some limitations and problems. It can only be applied to linear 
problems where the acoustic field beyond the boundary can be decomposed into 
an eigenvalue problem, and the eigenvalues and eigenmodes must be found at this 
location. In some geometries, they can be found analytically; however, for many 
geometries they must be found numerically. The boundary condition has been 
implemented in a panel method. It can also be applied to a finite element or 
finite difference method; however, this boundary condition will destroy the block 
tri-/penta- diagonal structure of the problem whenever more than one propagating 
mode is present. Thus, the solution procedure will have to be altered to allow a 
small, yet full, submatrix at any outflow boundary. 

SUMMARY 


Introducing artificial boundaries in acoustic computations can be difficult. The 
boundary condition must be very accurate in order to avoid reflections which are 
inconsistent with the actual acoustic field. To do this the boundary condition 
must represent the acoustic field outside the boundary. Simple boundary conditions 
specifying acoustic pressure, its normal derivative, or a combination of these will 
not work in most cases. 
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A boundary condition that matches the computation to a series representation 
of the outside acoustic field has been presented. A model problem shows that this 
boundary condition works well, even when several modes are present. 
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Table 1. - EIGENVALUES FOR TEST CASE WITH ONE PROPAGATING 
MODE, OUTFLOW BOUNDARY AT x = ±1.5. 


n m kj 


1 1 (l.,0.) 

1 2 (0.,4.36) 

2 1 (0.,2.97) 

2 2 (0.,5.37) 

3 1 (0.,6.19) 

3 2 (0.,7.64) 


A 

P 

analytic 

(0.7090,-0.0503) 

(-.00016,0.) 

(-.00154,0.) 

(-.000016,0.) 

(-.000018,0.) 

(-. 0000012 , 0 .) 


A 

P *i 

new b.c. 

(0.7043,-0.0835) 
(-. 0002 , . 0000 ) 
(-.0019, .0000) 
(. 0000 , . 0000 ) 
(.0015, .0050) 
(. 0000 , . 0000 ) 


A 

P *i 

simple b.c. 

(0.7036,-0.0720) 
(. 0000 , . 0000 ) 
(. 0000 , . 0000 ) 
(. 0000 ,. 0000 ) 
(. 0000 , . 0000 ) 
(. 0000 , . 0000 ) 
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Table 2. - EIGENVALUES FOR TEST CASE WITH THREE PROPAGATING 
MODES, OUTFLOW BOUNDARIES AT x = ±1.2. 


n 

m 

* i 

A 

Pd, 

analytic 

Pd, 

new b.c. 

A 

Pd, 

simple b.c. 

1 

1 

(6.200, 0.00000) 

(0.1049, -4.611E-02) 

(0.1048 -0.0486) 

(0.0431,-0.0927) 

1 

2 

(4.290, .00000) 

(-.1051, -.0491) 

(-0.1066 -0.0454) 

(.0000, .0000) 

1 

3 

(.0000, 6.455) 

(7.192E-05, .0000) 

(0.0008 0.0010) 

(.0000, .0000) 

1 

4 

(.0000, 11.90) 

(6.530E-08, .0000) 

(-0.0001 -0.0001) 

(.0000, .0000) 

2 

1 

(5.348, .0000) 

(9.988E-03, -7.332E-02) 

(0.0075 -0.0739) 

(.0000, .0000) 

2 

2 

(2.929, .00000) 

(-3.458E-02, 8.812E-02) 

(-0.0233 0.0873) 

(.0000, .0000) 

2 

3 

(.0000,7.176) 

(-1.514E-05, .0000) 

(0.0003 0.0000) 

(.oooo,.oooo) 

2 

4 

(.0000, 12.31) 

(-2.160E-08, .0000) 

(0.0000 0.0000) 

(.0000, .0000) 

3 

1 

(.00000, .9386) 

(-.4152, .0000) 

(-0.2927 -0.0738) 

(.0000, .0000) 

3 

2 

(.00000,4.572) 

(-7.620E-04, .0000) 

(-0.0008 -0.0004) 

(.0000, .0000) 

3 

3 

(.0000, 8.999) 

(-4.108E-06, .0000) 

(-0.0001 0.0009) 

(.0000, .0000) 

3 

4 

(.00000, 13.45) 

(-1.519E-08, .0000) 

(0.0000 0.0000) 

(.0000, .0000) 

4 

1 

(.0000, 7.073) 

(-3.098E-05, .0000) 

(0.0001 -0.0001) 

(.0000, .0000) 

4 

2 

(.0000, 8.370) 

(-3.872E-06, .0000) 

(0.0000 0.0001 ) 

(.0000, .0000) 

4 

3 

(.0000, 11.40) 

(-1.597E-07, .0000) 

(0.0000 0.0000) 

(.0000, .0000) 

4 

2 

(.0000, 15.17) 

(-1.523E-09, .0000) 

(0.0000 0.0000) 

(.0000, .0000) 

5 

1 

(.0000, 10.90) 

(-1.161E-07, .0000) 

(-0.0004 0.0002) 

(.0000, .0000) 

5 

2 

(.0000, 11.78) 

(-2.611E-08, .0000) 

(0.0000 0.0001) 

(.0000, .0000) 

5 

3 

(.0000, 14.10) 

(-2.900E-09, .0000) 

(0.0000 -0.0001) 

(.0000, .0000) 

5 

4 

(.0000, 17.29) 

(-6.000E-11, .0000) 

(.0000, .0000) 

(.0000, .0000) 
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FIGURES 



a) analytic solution, z = — .2 
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d)analytic solution, z = .2 




c) solution with simple boundary condi- 
tion, z = —.2 


f) solution with simple boundary condi- 
tion, z = .2 


Figure 2 - 9tt(p) for one propagating mode 
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Figure 3 — |p| for one propagating mode 
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Figure 4 - 5R(p) for five propagating modes 
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Figure 5 - |p| for five propagating modes 
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